#!/bin/bash 
#------------------------------------------------------------------------
# read ERA5 2d var data and draw 2D cape PIC
#set debug
 cmd=$0
 num_par=$#
#-----------------------------------------
echo "$cmd,$num_par"
if [ $num_par != 1 ] ;then
  echo " "
  echo "Error: no input parameter "
  echo " "
  exit
fi
varName=$1
# Directory information
#-----------------------------------------------------------------------
 HOMEDIR="//public/home/wangh29"
 MapDir="${HOMEDIR}/script/NCL_map/"
 DataDir="${HOMEDIR}/data2019-2021/script/paper/data"
 outdir="${HOMEDIR}/data2019-2021/script/paper/svg/newfig10ERA52dFRvsWR-${varName}"
 #----------
 #par
 start_lon=109. 
 end_lon=118. 
 start_lat=19.
 end_lat=26. 
# start_lon=70.
# end_lon=130.
# start_lat=10.
# end_lat=50.
 VINC=3
 cnLineThick=8
 FontS=0.025  # font size
 type="x11"  #pdf,x11,ncgm,ps
 RainID=14
 MinRainTh=25
 AWSrain="FALSE"
 RadarS="TRUE"
 MapLine="TRUE" #${MapLine}
 SecLine="FALSE"
#-----------------------------------------------------------------------
# Create plot.ncl 
#-----------------------------------------------------------------------
cat > comPressure.ncl << EOF
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl" 
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"

begin
;-----------------------------------------------------
; We generate plots, but what kind do we prefer?
     wks_type ="${type}"
   if("${type}".eq."png")then
    wks_type@wkWidth = 2500
    wks_type@wkHeight = 2500
   else
    wks_type@wkPaperHeightF  =11 
    wks_type@wkPaperWidthF = 11 
   end if 
   wks = gsn_open_wks(wks_type,"${outdir}")
;----------------------
;read ERA5
 if(1.eq.1)then
  FRfiles  = systemfunc(" csh -c | ls  "+ "${DataDir}/FR-FS.ERA5.2D/20*.nc")
  WRfiles  = systemfunc(" csh -c | ls  "+ "${DataDir}/WR-FS.ERA5.2D/20*.nc")
  ny = 201
  nx = 321
   gsn_define_colormap(wks,"MPL_Greens");"precip3_16lev")
 else
  FRfiles  = systemfunc(" csh -c | ls  "+ "${DataDir}/FR-FS.ERA5.2D/OtherVar/FR20*.nc")
  WRfiles  = systemfunc(" csh -c | ls  "+ "${DataDir}/WR-FS.ERA5.2D/OtherVar/WR20*.nc")
  ny = 101
  nx = 201
   gsn_define_colormap(wks,"GMT_no_green");"precip3_16lev")
 end if 

  nf = dimsizes(FRfiles)
  nw = dimsizes(WRfiles)
  nvar = 2
  cape = new((/nvar,ny,nx/),"double")
  cin = new((/nvar,ny,nx/),"double")
  cape = 0.
  cin  = 0.
  do ii=0,nf-1
   a    = addfile(FRfiles(ii),"r")
   time = wrf_user_getvar(a,"time",-1)
   lon  = wrf_user_getvar(a,"longitude",-1)     
   lat  = wrf_user_getvar(a,"latitude",-1)
   ce   = wrf_user_getvar(a,"${varName}",0)
   ci  = wrf_user_getvar(a,"tcslw",0)
   cape(0,:,:) = cape(0,:,:)+(ce*ce@scale_factor+ce@add_offset)/nf
   cin(0,:,:)  = cin(0,:,:)+(ci*ci@scale_factor+ci@add_offset)/nf
   delete([/ce,ci/])
  end do 
  binFile = "./data/FR-SurPress.bin"
  ;print(lon)
  ;print(lat)
  ;printVarSummary(cape)
setfileoption ("bin", "WriteByteOrder", "LittleEndian")
  fbinrecwrite(binFile,-1,lon)
  fbinrecwrite(binFile,-1,lat)
  fbinrecwrite(binFile,-1,cape(0,:,:))
  fbinrecwrite(binFile,-1,cin(0,:,:))

  do ii=0,nw-1
   b    = addfile(WRfiles(ii),"r")
   time = wrf_user_getvar(b,"time",-1)
   lon  = wrf_user_getvar(b,"longitude",-1)     
   lat  = wrf_user_getvar(b,"latitude",-1)
   ce   = wrf_user_getvar(b,"${varName}",0)
   ci  = wrf_user_getvar(b,"tcslw",0)
   cape(1,:,:) = cape(1,:,:)+(ce*ce@scale_factor+ce@add_offset)/nw
   cin(1,:,:)  = cin(1,:,:)+(ci*ci@scale_factor+ci@add_offset)/nw
   delete([/ce,ci/])
  end do
setfileoption ("bin", "WriteByteOrder", "LittleEndian")
  WRFile = "./data/WR-SurPress.bin"
  fbinrecwrite(WRFile,-1,lon)
  fbinrecwrite(WRFile,-1,lat)
  fbinrecwrite(WRFile,-1,cape(1,:,:))
  fbinrecwrite(WRFile,-1,cin(1,:,:))
;   cin = cin*1000. 
   lon@units      = "degrees-east"
   lat@units      = "degrees_north"
   cape!1 = "lattitude"
   cape!2 = "longttude"
   cape&longttude= lon
   cape&lattitude= lat
   cin!1 = "lattitude"
   cin!2 = "longttude"
   cin&longttude= lon
   cin&lattitude= lat
   tmpCape = cape(:,{${start_lat}:${end_lat}},{${start_lon}:${end_lon}})
   tmpCin  = cin(:,{${start_lat}:${end_lat}},{${start_lon}:${end_lon}})
   do i=0,nvar-1
    print("${varName}: "+min(tmpCape(i,:,:))+" "+max(tmpCape(i,:,:)))
    print("tcslw: "+min(tmpCin(i,:,:))+" "+max(tmpCin(i,:,:)))
   end do 

  res = True
;  res@vpWidthF = .65
;  res@vpHeightF= .7
  res@gsnAddCyclic= False
;  res@gsnMaximize = True      ; maximized in the workstation
  res@gsnDraw     = False
  res@gsnFrame    = False
  res@gsnLeftString=  "";"CAPE"; 
  res@gsnRightString = ""
  res@gsnCenterString  = "" 
;  res@gsnLeftStringOrthogonalPosF = -0.1
  res@gsnLeftStringFontHeightF   = ${FontS} 
  res@gsnRightStringFontHeightF  = ${FontS} 
  res@gsnCenterStringFontHeightF = ${FontS}
  res@mpOutlineOn         = True
  res@mpDataSetName       = "Earth..4"
  res@mpDataBaseVersion   = "MediumRes"
  res@mpOutlineSpecifiers = (/"Taiwan","Vietnam",\
                             "Thailand","Laos","Cambodia"/)
  res@mpFillColors        = (/"background","White","White","White", \
                              "transparent"/)
;Narrow the area to nedded Province
  res@mpLimitMode = "LatLon"
  res@mpMinLatF   =  ${start_lat}          ; nedded limits
  res@mpMaxLatF   =  ${end_lat}
  res@mpMinLonF   =  ${start_lon}
  res@mpMaxLonF   =  ${end_lon}


  res@mpGridAndLimbOn        = False 
  res@pmTickMarkDisplayMode  = "Always"
  res@mpGridSpacingF         = 10  ;
  res@mpGridLatSpacingF      = 10  ;
  res@mpGridLonSpacingF      =10  ;
  res@mpGridLineDashPattern  = 11   ;
 ;tickmark res
  res@tmXTOn                 = True;False
  res@tmYROn                 = True;False
  res@tmXTLabelsOn = False
  res@tmYRLabelsOn = False
  res@tmXBLabelsOn = True
  res@tmYLLabelsOn = True
  res@tmXBLabelStride = 2
  res@tmYLLabelStride = 2
  res@tmXBLabelFontHeightF  = ${FontS} ;0.015
  res@tmYLLabelFontHeightF  = ${FontS} ;0.015


;--------------------------------
; plot for divergence
;--------------------------------
  opts_d = res
  opts_d@cnFillOn     = True 
  opts_d@cnLinesOn    = False
  opts_d@cnLineLabelsOn       = False ; no contour line labels
  opts_d@cnInfoLabelOn= False
  opts_d@cnLevelSelectionMode  =  "Explicitlevels"
  Titles = ""
  if("${varName}".eq."cape")then
   opts_d@cnLevels  = ispan(500,3500,250)
   Titles = "CAPE (J kg ~S~-1~N~)"
  elseif("${varName}".eq."cin")then
   opts_d@cnLevels  = ispan(0,600,50)
   Titles = "CIN (J kg ~S~-1~N~)"
  elseif("${varName}".eq."tcwv")then
   opts_d@cnLevels  = ispan(44,64,2)
   Titles = "Column vapour (kg m ~S~-2~N~)"
  elseif("${varName}".eq."tcw")then
   opts_d@cnLevels  = ispan(44,64,2)
   Titles = "Column water (kg m ~S~-2~N~)"
  elseif("${varName}".eq."tcrw")then
   opts_d@cnLevels  = ispan(44,64,2)
  elseif("${varName}".eq."tcslw")then
   opts_d@cnLevels  = fspan(0.,0.15,16)
   Titles = "Column supercooled liquid water(kg m ~S~-2~N~)"
  elseif("${varName}".eq."totalx")then
   opts_d@cnLevels  = ispan(38,46,1)
   Titles = "totals index (K)"
  elseif("${varName}".eq."kx")then
   opts_d@cnLevels  = ispan(22,38,1)
   Titles = "K index (K)"
  elseif("${varName}".eq."cp")then
   cape = cape*1000.
   opts_d@cnLevels  = fspan(0.,4.,11)
   Titles = "Convective precipitation (mm)"
  elseif("${varName}".eq."mvimd")then
   cape = cape*10000.
   opts_d@cnLevels  = ispan(-10,10,2);fspan(0.,4.,11)
   Titles = "Integrated moisture divergence (10~S~-4~N~ kg m~S~-2~N~ s~S~-1~N~)"
  elseif("${varName}".eq."vimd")then
   opts_d@cnLevels  = fspan(-4.,4.,9);fspan(-3.,3.,13)
   Titles = "Vertically integrated moisture divergence ( kg m~S~-2~N~)"
  elseif("${varName}".eq."sp")then
   cape = cape/100.
   opts_d@cnLevels  = ispan(894,1008,4);fspan(-3.,3.,13)
   Titles = "Surface pressure (hPa)"
  end if 

  opts_d@lbLabelBarOn       = False
  opts_d@lbAutoManage       = False             ; we control label bar
  opts_d@pmLabelBarWidthF   = 0.62
  opts_d@pmLabelBarHeightF     = 0.1               ; default is taller
  opts_d@lbLabelFontHeightF    = ${FontS}*0.7;         ; default is HUGE
  opts_d@lbPerimOn             = False             ; default has box
  opts_d@lbTopMarginF       = 0.3

  plotA = new(nvar,graphic)
  plotA(0) = gsn_csm_contour_map(wks,cape(0,:,:),opts_d)
  opts_d@tmYLLabelsOn = False 
  plotA(1) = gsn_csm_contour_map(wks,cape(1,:,:),opts_d)

;---------------------------------
; Plot for hgt 
;---------------------------------
    opts_t = True
    opts_t@gsnDraw      =  False
    opts_t@gsnFrame     =  False
    opts_t@gsnLeftString  =  ""
    opts_t@gsnRightString =  ""
    opts_t@cnFillOn	= False 
    opts_t@cnLinesOn	= True
    opts_t@cnLineColor	= "red"
    opts_t@cnInfoLabelOn= False
    opts_t@cnLevelSelectionMode  =  "Explicitlevels"
    opts_t@cnLevels  =  fspan(0.04,0.14,3) ;ispan(0,600,100)
; For Line Label
    opts_t@cnLabelMasking             = True
    opts_t@cnLineLabelFontColor       = opts_t@cnLineColor
    opts_t@cnLineLabelBackgroundColor = "transparent"
;   opts_t@cnLineLabelPlacementMode   = "Constant"
;   opts_t@cnLineDrawOrder           = "PostDraw"
    opts_t@cnLineLabelFontHeightF     = ${FontS}*0.8;.012
    opts_t@cnInfoLabelOrthogonalPosF  = 0.1
    opts_t@cnLineLabelInterval        = 1 
    opts_t@cnLineLabelDensityF       = 1   ;
    opts_t@cnInfoLabelOrthogonalPosF  = 0.3  ; offset second label information
    opts_t@cnLineLabelPerimOn = False
    opts_t@cnInfoLabelOn     = False

  ;  opts_t@cnInfoLabelOrthogonalPosF = 0.07  ; offset second label information
    opts_t@gsnContourLineThicknessesScale = 2.0
    plotD = new(nvar,graphic)
    do i=0,nvar-1
      cin(i,:,:) = smth9(cin(i,:,:),0.5,-0.25,False)
      plotD(i) = gsn_csm_contour(wks,cin(i,:,:),opts_t)
      overlay(plotA(i),plotD(i))

   end do 
;---------------------------------------------------
  if("${MapLine}".eq."TRUE")then
   mapdir = "${MapDir}"
   shp_city=mapdir+"cantonnew2.shp"
   shp_province=mapdir+"gdnew.shp"
   shp_huanan=mapdir+"SC_province.shp"
   shp_state=mapdir+"CHINA.shp"
   pres               = True
   pres@cnFillOn      = False
   pres@gsLineThicknessF=2
   pres@gsLineColor   = "grey33"
   poly_2 = new((/nvar,51/),graphic) 
   poly_3 = new((/nvar,1288/),graphic)
   do i=0,nvar-1
  ;  poly_2(i,:)=gsn_add_shapefile_polylines(wks, plotA(i), shp_province, pres)
    poly_3(i,:)=gsn_add_shapefile_polylines(wks, plotA(i), shp_state, pres)
   end do 
  end if 
;-----------------------------------------
; radar site 
 if("${RadarS}".eq."TRUE")then 
  dpr_lon = (/113.37,111.9789/)
  dpr_lat = (/23.0039,21.845/)
  pmres = True
  pmres@gsMarkerSizeF = 0.015
  pmres@gsMarkerOpacityF = 1.0
  pmres@gsMarkerColor = "dim gray"
  pmres@gsMarkerIndex = 12 
  pmres@gsMarkerThicknessF = 5
  radarM = new(nvar,"graphic")
  ;-------------------------
  ;Add radar 150 km OBS range
  pi    = 3.1415926
  an2ra = pi/180.
  radius= 150./(6371*an2ra) ;2.069
  nRad  = dimsizes(dpr_lon)
  poly = new((/nvar,nRad/),"graphic")
  X_cir= new((/nRad,400/),"float")
  Y_cir= new((/nRad,400/),"float")
  theta = fspan(0,6.28,400)
  c     = cos(theta)
  s     = sin(theta)
  do j=0,nRad-1
   do i = 0 ,400-1
      X_cir(j,i) = dpr_lon(j) + radius*c(i)/cos(dpr_lat(j)*an2ra)
      Y_cir(j,i) = dpr_lat(j) + radius*s(i)
   end do
  end do
  pl    = True
  pl@gsLineThicknessF = 4            ; line thickness
  pl@gsLineColor      = pmres@gsMarkerColor 
  pl@gsLineDashPattern= 2
  do i=0,nvar-1
   radarM(i) = gsn_add_polymarker(wks,plotA(i),dpr_lon,dpr_lat,pmres) 
   do j=0,nRad-1
     poly(i,j) = gsn_add_polyline(wks,plotA(i),X_cir(j,:),Y_cir(j,:),pl)
   end do
  end do 
 end if 
;----------------------------------------------
   pnlres                         = True
;   pnlres@pmTitleDisplayMode      = "Always" 
   pnlres@txString                = ""
   pnlres@gsnPanelFigureStrings   = (/"a) FR","b) WR"/)
   pnlres@gsnPanelFigureStringsFontHeightF = ${FontS}*0.7 
   pnlres@amJust = "BottomRight";"TopLeft"
   pnlres@gsnPanelXF  = (/0.11,0.53/)  ; slightly adjust X location of each plot

   pnlres@gsnPanelFigureStringsPerimOn = True
   pnlres@gsnPanelLabelBar           = True    ; Turn on common labelbar
   pnlres@lbLabelAutoStride          = True    ; Spacing of lbar labels.
   pnlres@lbTopMarginF        = 0.36
   pnlres@lbTitleString       = Titles
   pnlres@lbTitlePosition    = "Bottom"
   pnlres@lbTitleDirection   = "Across"
   pnlres@lbTitleFontHeightF = ${FontS}*0.8
   pnlres@lbTitleOffsetF     = 0.06
   pnlres@pmLabelBarWidthF   = 0.7
   pnlres@pmLabelBarHeightF  = 0.08               ; default is taller
   pnlres@lbLabelFontHeightF    = ${FontS}*0.8  ; default is HUGE

   pnlres@gsnMaximize     = True               ; 
;   pnlres@gsnPaperOrientation = "auto"
;   pnlres@gsnPaperWidth   = 15. 
   gsn_panel(wks,plotA,(/1,2/),pnlres)


end
EOF

ncl comPressure.ncl
rm -f comPressure.ncl 
exit
